﻿"""
Jieun Oh

mdct.py -- Jieun's implementation of MDCT for Music 422 HW3 #1

"""
import numpy as np


#------------------------------MDCT slow-----------------------------#
def MDCTslow(data, a, b, isInverse=False):
     """slow MDCT & IMDCT implementation using p.130 """
     N = a+b
     n0 = (b+1)/2.

     r = N
     c = N/2
    
     matrix = np.zeros((r,c))

     for i in range(r):
        for j in range(c):
           matrix[i,j] = np.cos(2*np.pi*(i+n0)*(j+0.5)/N)
 
     if (isInverse):
        matrix = (matrix.T)   

     else:
        matrix = matrix* (4./N)    # multiply by 4/N for forward transform only

     return np.dot(data, matrix)



#------------------------------FFT-based MDCT----------------------------#

def MDCT(data, a, b, isInverse=False):
     """fast fft-based MDCT & IMDCT implementation using p.141-143 """

     N = a+b              #np.size(data)
     n0 = (b+1)/2.

     if (not isInverse):

        # pre-twiddle
        nIndices = np.arange(N)
        preMultiplyArray = np.cos(np.pi*nIndices/N) + np.sin(-np.pi*nIndices/N)*1j
        data = data*preMultiplyArray

        # N-point FFT on pre-twiddled data
        #####transform = np.fft.fft(data)                ## try Isaac's approach: *2.0/N
        transform = np.fft.fft(data) * 2.0/N                ## try Isaac's approach: *2.0/N

        # post-twiddle
        kIndices = np.arange(N/2)
        post = transform[0:N/2]
        postMultiplyArray = np.cos(2*np.pi*n0*(kIndices+0.5)/N) + np.sin(-2*np.pi*n0*(kIndices+0.5)/N)*1j
        #####answer = np.real(post*postMultiplyArray)*2  ## try Isaac's approach: No *2
        answer = np.real(post*postMultiplyArray)  ## try Isaac's approach: No *2

     else:

        # pre-twiddle
        negativeFlip = -1*data[::-1]
        data = np.concatenate((data, negativeFlip))
        kIndices = np.arange(N)
        preMultiplyArray = np.cos(2*np.pi*n0*kIndices/N) + np.sin(2*np.pi*n0*kIndices/N)*1j
        data = data*preMultiplyArray

        # N-point IFFT on pre-twiddled data
        #####transform = np.fft.ifft(data)              ## try Isaac's approach: *N
        transform = np.fft.ifft(data) *N             ## try Isaac's approach: *N

        # post-twiddle
        nIndices = np.arange(N)
        post = transform
        postMultiplyArray = np.cos(np.pi*(nIndices+n0)/N) + np.sin(np.pi*(nIndices+n0)/N)*1j
        answer = np.real(post*postMultiplyArray) 


     return answer





def IMDCT(data, a, b):
     return MDCT(data, a, b, True)



# --------------------------Testing the MDCT--------------------------------#


if __name__=="__main__":

   from timeit import Timer

   def slowMDCT():
 
      b1 = np.array([0,0,0,0,0,1,2,3])/np.sqrt(2)
      b2 = np.array([0,1,2,3,4,5,6,7])/np.sqrt(2)
      b3 = np.array([4,5,6,7,8,9,10,11])/np.sqrt(2)
      b4 = np.array([8,9,10,11,0,0,0,0])/np.sqrt(2)

      a = 4
      b = 4
      
      '''
      # uncomment this section to test execution timing on a long array
      b1 = np.concatenate((np.zeros(250), np.arange(0, 250)))
      b2 = np.arange(0, 500)
      b3 = np.arange(250, 750)
      b4 = np.concatenate((np.arange(750, 1000), np.zeros(250)))
      
      a = 250
      b = 250
      '''

      fwd_b1 = MDCTslow( b1, a, b)
      fwd_b2 = MDCTslow( b2, a, b)
      fwd_b3 = MDCTslow( b3, a, b)
      fwd_b4 = MDCTslow( b4, a, b)

      inv_b1 = MDCTslow( fwd_b1, a, b, True)/np.sqrt(2)
      inv_b2 = MDCTslow( fwd_b2, a, b, True)/np.sqrt(2)
      inv_b3 = MDCTslow( fwd_b3, a, b, True)/np.sqrt(2)
      inv_b4 = MDCTslow( fwd_b4, a, b, True)/np.sqrt(2)

      zeros = np.zeros(a)
      new_b = np.zeros(3*a)
      new_b[0:a] = inv_b1[a:2*a] + inv_b2[0:a]
      new_b[a:2*a] = inv_b2[a:2*a] + inv_b3[0:a]
      new_b[2*a:3*a] = inv_b3[a:2*a] + inv_b4[0:a]

      return new_b


   def fastMDCT():

      b1 = np.array([0,0,0,0,0,1,2,3])/np.sqrt(2)
      b2 = np.array([0,1,2,3,4,5,6,7])/np.sqrt(2)
      b3 = np.array([4,5,6,7,8,9,10,11])/np.sqrt(2)
      b4 = np.array([8,9,10,11,0,0,0,0])/np.sqrt(2)

      a = 4
      b = 4

      '''
      # uncomment this section to test execution timing on a long array
      b1 = np.concatenate((np.zeros(250), np.arange(0, 250)))
      b2 = np.arange(0, 500)
      b3 = np.arange(250, 750)
      b4 = np.concatenate((np.arange(750, 1000), np.zeros(250)))
      
      a = 250
      b = 250
      '''

      fwd_b1 = MDCT( b1, a, b)
      fwd_b2 = MDCT( b2, a, b)
      fwd_b3 = MDCT( b3, a, b)
      fwd_b4 = MDCT( b4, a, b)
   
      inv_b1 = MDCT( fwd_b1, a, b, True)/np.sqrt(2)
      inv_b2 = MDCT( fwd_b2, a, b, True)/np.sqrt(2)
      inv_b3 = MDCT( fwd_b3, a, b, True)/np.sqrt(2)
      inv_b4 = MDCT( fwd_b4, a, b, True)/np.sqrt(2)

      zeros = np.zeros(a)
      new_b = np.zeros(3*a)
      new_b[0:a] = inv_b1[a:2*a] + inv_b2[0:a]
      new_b[a:2*a] = inv_b2[a:2*a] + inv_b3[0:a]
      new_b[2*a:3*a] = inv_b3[a:2*a] + inv_b4[0:a]
  
      return new_b


   t = Timer("slowMDCT()", "from __main__ import slowMDCT")
   print t.timeit(number=100)
   # results in approx. 13.699s for length 12 array, repeating 10000 times
   # results in approx. 42.416s for length 750 array, repeating 10 times

   t = Timer("fastMDCT()", "from __main__ import fastMDCT")
   print t.timeit(number=100)
   # results in approx. 10.480s for length 12 array, repeating 10000 times
   # results in approx. 0.0212s for length 750 array, repeating 10 times


   # new_b = slowMDCT()
   # new_b = fastMDCT()
   # print(new_b)




""" Answers to #1 """

# Description of aliased portion of the IMDCT output:
#    --> inv_b1 = [0, 0, 0, 0, 1.5, 1.5, 1.5, 1.5]
#        inv_b2 = [-1.5, -0.5, 0.5, 1.5, 5.5, 5.5, 5.5, 5.5]
#        inv_b3 = [-1.5, -0.5, 0.5, 1.5, 9.5, 9.5, 9.5, 9.5]
#        inv_b4 = [-1.5, -5, 5, 1.5, 0, 0, 0, 0]
#    --> The result of MDCT followed by IMDCT results in severe aliasing; 
#        The output is quite different from the input. 
# Describe how OLA can exactly recover the input signal:
#    --> However, overlap-adding the results perfectly recovers the signal because
#        the two halves make-up for each other's deficiencies. 
# Flat window of 1/sqrt(2) satisfy the perfect reconstruction window condition: 
#    --> (doubling from 50% overlap) * analysis window * synthesis window
#         = 2*(1/sqrt(2))*(1/sqrt(2)) = 1. 


# Comparing Execution Time of slow vs. fast implementation of MDCT/ IMDCT:
#    --> The speed becomes very noticeable for longer arrays. 
#        As noted above, for length 750 arrays, the slow implementation repeated 10 times takes 
#        approx. 42.416s, in comparison to 0.0212s for the faster implementation.  



